{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OpenSky flight trajectories\n",
    "\n",
    "Flight path information for commercial flights is available for some regions of the USA and Europe from the crowd-sourced [OpenSky Network](https://opensky-network.org/).  OpenSky collects data from a large number of users monitoring public air-traffic control information.  Here we will use a subset of the data that was polled from their REST API at an interval of 1 minute over 4 days (September 5-13, 2016), using the scripts shown at the end of this notebook.  In general the terms of use for OpenSky data do not allow redistribution, but we have obtained specific permission for distributing a[this subset of the data](https://s3.amazonaws.com/datashader-data/opensky.parq), which is a 200MB Parquet file (1.1GB as the original database). If you want more or different data, you can run the scripts at the end of this notebook to collect some yourself, or else you can contact OpenSky asking for a copy of the dataset.\n",
    "\n",
    "We'll only use some of the fields provided by OpenSky, out of: *icao24, callsign, origin, time_position, time_velocity, longitude, latitude, altitude, on_ground, velocity, heading, vertical_rate, sensors, timestamp*\n",
    "\n",
    "Here, we'll load the data and declare that some fields are categorical (which isn't information fully expressed in the Parquet file):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "import pandas as pd\n",
    "\n",
    "flightpaths = pd.read_parquet('../data/opensky.parq')\n",
    "flightpaths['origin']    = flightpaths.origin.astype('category')\n",
    "flightpaths['ascending'] = flightpaths.ascending.astype('category')\n",
    "flightpaths.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The default database has about 10 million points, with some metadata for each.  \n",
    "\n",
    "Now let's define a datashader-based processing pipeline to render images:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datashader as ds\n",
    "import datashader.transfer_functions as tf\n",
    "from colorcet import fire\n",
    "from matplotlib.colors import rgb2hex\n",
    "from matplotlib.cm import get_cmap\n",
    "\n",
    "import numpy as np\n",
    "from cartopy import crs\n",
    "\n",
    "plot_width  = 850\n",
    "plot_height = 600\n",
    "x_range = (-2.0e6, 2.5e6)\n",
    "y_range = (4.1e6, 7.8e6)\n",
    "\n",
    "def categorical_color_key(ncats,cmap):\n",
    "    \"\"\"Generate a color key from the given colormap with the requested number of colors\"\"\"\n",
    "    mapper = get_cmap(cmap)\n",
    "    return [str(rgb2hex(mapper(i))) for i in np.linspace(0, 1, ncats)]\n",
    "\n",
    "def create_image(x_range=x_range, y_range=y_range, w=plot_width, h=plot_height, \n",
    "                 aggregator=ds.count(), categorical=None, black=False, cmap=\"blue\"):\n",
    "    opts={}\n",
    "    if categorical and cmap:\n",
    "        opts['color_key'] = categorical_color_key(len(flightpaths[aggregator.column].unique()),cmap)       \n",
    "\n",
    "    cvs = ds.Canvas(plot_width=w, plot_height=h, x_range=x_range, y_range=y_range)\n",
    "    agg = cvs.line(flightpaths, 'longitude', 'latitude',  aggregator)\n",
    "    img = tf.shade(agg, cmap=cmap, **opts)\n",
    "        \n",
    "    if black: img = tf.set_background(img, 'black')\n",
    "    return img"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use this function to get a dump of all of the trajectory information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "create_image(aggregator=ds.count(), cmap=fire, black=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This plot shows all of the trajectories in this database, overlaid in a way that avoids [overplotting](https://anaconda.org/jbednar/plotting_pitfalls/notebook).  With this \"fire\" color map, a single trajectory shows up as black, while increasing levels of overlap show up as brighter colors.  \n",
    "\n",
    "A static image on its own like this is difficult to interpret, but if we overlay it on a map we can see where these flights originate, and can zoom in to see detail in specific regions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader.bokeh_ext import InteractiveImage\n",
    "from bokeh.plotting import figure, output_notebook\n",
    "from bokeh.models.tiles import WMTSTileSource\n",
    "\n",
    "output_notebook()\n",
    "\n",
    "def base_plot(tools='pan,wheel_zoom,reset',plot_width=plot_width, plot_height=plot_height,**plot_args):\n",
    "    p = figure(tools=tools, plot_width=plot_width, plot_height=plot_height,\n",
    "        x_range=x_range, y_range=y_range, outline_line_color=None,\n",
    "        min_border=0, min_border_left=0, min_border_right=0,\n",
    "        min_border_top=0, min_border_bottom=0, **plot_args)\n",
    "    \n",
    "    p.axis.visible = False\n",
    "    p.xgrid.grid_line_color = None\n",
    "    p.ygrid.grid_line_color = None\n",
    "    \n",
    "    return p\n",
    "\n",
    "ArcGIS=WMTSTileSource(url='http://server.arcgisonline.com/ArcGIS/rest/services/'\n",
    "                      'World_Street_Map/MapServer/tile/{Z}/{Y}/{X}.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = base_plot()\n",
    "p.add_tile(ArcGIS)\n",
    "InteractiveImage(p, create_image, aggregator=ds.count())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "E.g. try zooming in on London in the above figure, which has a lot of structure not visible in the initial rendering but visible on a zoom. Note that zooming in will only reveal more detail in the datashader plot if you are working with a live server; a static HTML view (e.g. on Anaconda Cloud) will dynamically update the underlying map plot, but not the data.  \n",
    "\n",
    "We can use the metadata associated with each trajectory to show additional information.  For instance, we can color each flight by its country of origin, using the key:\n",
    "\n",
    "* **UK** - Orange\n",
    "* **Germany** - Blue\n",
    "* **Netherland** - Teal\n",
    "* **Switzerland** - Yellow\n",
    "* **France** - Purple\n",
    "* **Norway** - Green\n",
    "* **USA** - Red\n",
    "\n",
    "(There are actually more than a hundred different origins, so this key is only approximate.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = base_plot()\n",
    "p.add_tile(ArcGIS)\n",
    "InteractiveImage(p, create_image, categorical=True, aggregator=ds.count_cat('origin'), cmap='hsv_r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can label ascending (Blue) vs. descending flights (Red), which is particularly informative when zooming in on specific airports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = base_plot()\n",
    "p.add_tile(ArcGIS)\n",
    "InteractiveImage(p, create_image, aggregator=ds.count_cat('ascending'), cmap=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can show velocity, which of course decreases (dark colors) when approaching or leaving airports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = base_plot()\n",
    "p.add_tile(ArcGIS)\n",
    "InteractiveImage(p, create_image, aggregator=ds.mean('velocity'), cmap=fire[::-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The flight patterns associated with each airport are clearly visible in these close-ups of various cities, where the circular holding pattern for landings (red) is clearly visible for the various airports in London:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import holoviews as hv\n",
    "from holoviews import opts\n",
    "hv.extension('matplotlib')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=140)\n",
    "opts.defaults(\n",
    "    opts.RGB(xaxis=None, yaxis=None),\n",
    "    opts.Layout(hspace=0.1, vspace=0, sublabel_format=None)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_rgb(img):\n",
    "    return np.flipud(img.view(dtype=np.uint8).reshape(img.shape[:2] + (4,)))\n",
    "\n",
    "cities = {'Frankfurt' : (8.6821, 50.1109),\n",
    "          'London'    : (-0.1278, 51.5074), \n",
    "          'Paris'     : (2.3522, 48.8566),\n",
    "          'Amsterdam' : (4.8952, 52.3702),\n",
    "          'Zurich'    : (8.5417, 47.3769),\n",
    "          'Munich'    : (11.5820, 48.1351)}\n",
    "\n",
    "radius = 150000\n",
    "\n",
    "mercator_cities = {city: crs.GOOGLE_MERCATOR.transform_point(lon, lat, crs.PlateCarree()) \n",
    "                   for city, (lon, lat) in cities.items()}\n",
    "city_ranges = {city: dict(x_range=(lon-radius, lon+radius), y_range=(lat-radius, lat+radius))\n",
    "               for city, (lon, lat) in mercator_cities.items()}\n",
    "\n",
    "hv.Layout([hv.RGB(to_rgb(create_image(aggregator=ds.count_cat('ascending'), black=True, categorical=True, \n",
    "                                      w=300, h=300, cmap=None, **ranges).data), group=city)\n",
    "                    for city, ranges in sorted(city_ranges.items())]).cols(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or colorized by flight origin:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.Layout([hv.RGB(to_rgb(create_image(aggregator=ds.count_cat('origin'), black=True,\n",
    "                                      categorical=True, w=300, h=300, cmap='hsv_r', **ranges).data), group=city)\n",
    "                    for city, ranges in sorted(city_ranges.items())]).cols(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The patterns for a single city can make a nice wallpaper for your desktop if you wish:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "city,ranges = \"Zurich\",city_ranges[\"Zurich\"]\n",
    "create_image(aggregator=ds.count_cat('origin'), black=False,\n",
    "                           categorical=True, w=800, h=800, cmap='hsv_r', **ranges)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, datashader makes it quite easy to explore even large databases of trajectory information, without trial and error parameter setting and experimentation.  These examples have millions of datapoints, but it could work with [billions](http://anaconda.org/jbednar/osm/notebook) just as easily, covering long time ranges or large geographic areas. Check out the other [datashader notebooks](http://anaconda.org/jbednar/notebooks) for other examples!\n",
    "\n",
    "\n",
    "## Downloading and preparing the data\n",
    "\n",
    "This data was obtained by running a cron job with the following script running at one-minute intervals over a four-day period:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "import json, sqlite3, requests, pandas as pd\n",
    "\n",
    "DB='../data/opensky.db'\n",
    "conn = sqlite3.connect(DB)\n",
    "api_url = 'https://opensky-network.org/api/states/all'\n",
    "\n",
    "cols = ['icao24', 'callsign', 'origin', 'time_position', 'time_velocity', \n",
    "        'longitude', 'latitude', 'altitude', 'on_ground', 'velocity', \n",
    "        'heading', 'vertical_rate', 'sensors']\n",
    "\n",
    "req = requests.get(api_url)\n",
    "content = json.loads(req.content)\n",
    "states = content['states']\n",
    "df = pd.DataFrame(states, columns=cols)\n",
    "df['timestamp'] = content['time']\n",
    "df.to_sql('opensky', conn, index=False, if_exists='append')\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting `opensky.db` file was then transformed into Web Mercator coordinates, split per flight, and exported to Parquet format, using the code below.  This process took about 7 minutes on a MacBook Pro laptop."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "import sqlite3, pandas as pd, numpy as np, holoviews as hv, datashader.utils as du\n",
    "\n",
    "def transform_coords(df):\n",
    "    df=df.copy()\n",
    "    df.loc[:, 'longitude'], df.loc[:, 'latitude'] = \\\n",
    "        du.lnglat_to_meters(df.longitude,df.latitude)\n",
    "    return df\n",
    "\n",
    "def split_flights(df):\n",
    "    df = df.copy().reset_index(drop=True)\n",
    "    df = df[np.logical_not(df.time_position.isnull())]\n",
    "    empty=df[:1].copy()\n",
    "    empty.loc[0, :] = 0\n",
    "    empty.loc[0, 'origin'] = ''\n",
    "    empty.loc[0, 'latitude'] = np.NaN\n",
    "    empty.loc[0, 'longitude'] = np.NaN\n",
    "    paths = []\n",
    "    for gid, group in df.groupby('icao24'):\n",
    "        times = group.time_position\n",
    "        splits = np.split(group.reset_index(drop=True), np.where(times.diff()>600)[0])\n",
    "        for split_df in splits:\n",
    "            if len(split_df) > 20:\n",
    "                paths += [split_df, empty]\n",
    "    split = pd.concat(paths,ignore_index=True)\n",
    "    split['ascending'] = split.vertical_rate>0\n",
    "    return split\n",
    "\n",
    "# Load the data from a SQLite database and project into Web Mercator (1.5 min)\n",
    "DB='../data/opensky.db'\n",
    "conn = sqlite3.connect(DB)\n",
    "df = transform_coords(pd.read_sql(\"SELECT * from flights\", conn))\n",
    "\n",
    "# Split into groups by flight (6 min)\n",
    "flightpaths = split_flights(df)\n",
    "\n",
    "# Remove unused columns and declare categoricals\n",
    "flightpaths = flightpaths[['longitude', 'latitude', 'origin', 'ascending', 'velocity']]\n",
    "flightpaths['origin']    = flightpaths.origin.astype('category')\n",
    "flightpaths['ascending'] = flightpaths.ascending.astype('bool')\n",
    "\n",
    "# Export to Parquet\n",
    "args = dict(engine=\"fastparquet\", compression=\"snappy\", has_nulls=False, write_index=False)\n",
    "flightpaths.to_parquet(\"../data/opensky.parq\", **args)\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
